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Abstract —This paper presents a hybrid differential evolution 
with multiplier updating (HDE-MU) to solve the complex 
combined heat and power economic dispatch (CHPED) 
problems. Transmission losses and valve-point effects of 
conventional thermal generators are taken into account. The 
hybrid differential evolution (HDE) has the ability to efficiently 
search and actively explore solutions. Multiplier updating (MU) 
is introduced to avoid deforming the augmented Lagrange 
function (ALF), which is adopted to manage system constraints 
of the CHPED problem. The proposed HDE-MU integrates the 
HDE with the MU. A practical CHPED system is employed to 
demonstrate that the proposed algorithm has the benefits of 
straightforwardness; ease of implementation; better 
effectiveness than the previous methods, and the requirement 
for only a small population when applied to the CHPED 
operation. 

Index Terms —Combined heat and power, differential 
evolution, economic dispatch, multiplier updating. 


I. INTRODUCTION 

The conversion of fossil fuel into electricity is not an 
efficient process. Even the energy efficiency of the most 
modern combined cycle plants are between 50 and 60% 
efficient [1], Most of the energy wasted in the conversion 
process is heat. Recently, combined heat and power (CHP) 
units have played an increasingly important role in the utility 
industry [2~ 6]. Complexity arises if one or more units 
produce both power and heat. Such a case, both the power 
demand and the heat demand must be satisfied. Cogeneration 
units provide both electrical power and heat to customers. The 
heat production capacity of most cogeneration units, depend 
on the power generated, and vice versa. The mutual 
dependency of multiple-demand and Heat-Power capacity of 
those units introduce complexities in the integration of 
cogeneration units into the economic dispatch problem. The 
CHPED problem even may be more complicated if 
transmission losses and valve-point effects are taken into 
account. Solving such a complex optimization problem 
requires powerful techniques. Non-linear optimal algorithms, 
such as dual and quadratic programming [7], and gradient 
descent approaches, such as Lagrangian relaxation [8], have 
been applied for solving CHPED problems. However, these 
methods cannot handle non-smooth non-convex fuel cost 
function of the conventional thermal generator. The advent of 
stochastic search algorithms has overcome this problem for 
solving CHPED problems. Gravitational search algorithm 
(GSA) [2], cuckoo search algorithm [3], group search 
optimization (GSO) [4], exchange market algorithm [5], real 
coded genetic algorithm [6], opposition-based group search 
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optimization (OGSO) [9], Improved ant colony search 
algorithm [10], evolutionary programming [11], genetic 
algorithm [12], harmonic search algorithm [13], hybrid 
particle swarm optimization [14], self-adaptive real-coded 
genetic algorithm [15], novel selective particle swarm 
optimization [16], mesh adaptive direct search algorithm [17], 
oppositional teaching learning based optimization (OTLBO) 
[18], particle swarm optimization with time varying 
acceleration coefficients (TVAC-PSO) [19] and krill herd 
algorithm (KHA) [20] have been applied to solve the CHPED 
problems. The results obtained by the proposed algorithm, 
involving HDE and MU, have been compared with the 
existing methodologies reported in recent literature. 

Storn and Price [21] developed the Differential evolution 
(DE) which immediately gained popularity as a robust 
evolutionary algorithm. DE has been widely applied to the 
optimization problems in the power systems [22~ 26]. 
Throughout the years, DE has been used extensively for 
optimization problems, many results of which are the best 
compared to other standard methodologies. However, it has 
problems of converging onto local optimal solutions. A 
migration [27] is embedded in the proposed HDE to 
overcome such drawbacks for solving CHPED problems. The 
migration operation is included in the HDE to regenerate a 
newly diverse population, which prevents individuals from 
gradually clustering and thus greatly increases the amount of 
search space explored for a small population. 

The rest of this paper is organized as follows: in the next 
section, the basic equations governing CHPED have been 
briefly outlined. Section III explains the basic framework of 
the proposed algorithm involving HDE and MU. Section IV 
describes the test system, compares the results obtained by 
HDE-MU with other methods and provides the analysis of the 
obtained results. Finally, the conclusion is drawn in Section 
V. 

II. Problem formulation 
A. Objective function 

The CHPED problem is to determine the unit power and 
heat production so that the system’s production cost is 
minimized while the power and heat demands and other 
constraints are met. It can be mathematically stated as, 

n p n b , , 

Minimize ^ Cost, (/])+1 Cost j [H j , P,)+ X Cost, (H k ) (1) 

1=1 7=1 k =1 

Where 

Cost i (f) = a I + b,P, + c,P, 2 + |e ; x sin(/, x (/> min - p. | (2) 

Costj(Hj Pj) = cij +bjPj +CjPf +djH 2 +>'jH i +WjPjHj( 3) 
Cost k (H k ) = a k + b k H k + c k H k (4) 
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and where P stands unit power generation; H is unit heat 
production; Cost HP ^ with i=l,2,..., n p represents the cost 
function of the ; th power-only unit in $/h and a„ /;„ and c,- are 
cost coefficients of generator i, and e, and f are fuel cost 
coefficients of unit i with valve-points effects; Costj(Hj,Pj) 
withy=7,2,.... n b is the fuel cost function of the/ h CHP unit in 
$/h and a p bj, Cj, dj, rj, and wj are cost coefficients of unit j; 
CostfHf) with k=l,2,..., n h represents the cost function of the 
k th heat-only unit in $/h and a k , b k . and c k are cost coefficients; 
n p , n b , and n h denote the number of power-only units, 
co-generation units, and heat-only units, respectively; 
subscripts of i, j, and k are used for above mentioned unit; it 
should be mentioned that (P„ Pj) and (//,, H k ) are the output 
active power (in MW) and heat production (in MWth), 
respectively. 

B. Constraints 

Subject to the equilibrium constraints of electricity and 
heat production, and the capacity limits of each unit. Power 
and heat balance equations are represented by equality 
constraints (5) and (6), respectively. 


N p vectors of variables x in the optimization problem, namely, 
x G = \x G ,i = l,...,/V /( j, as a population in generation G. For 

convenience, the decision vector, x h is represented as 
(x u ...x l ,...x„ ,)■ Here, the decision variable, x Jt is directly 

coded as a real value within its bounds of (x™ n , x” 3 *). Each 
individual is generated as follows: 

x?| =xf n +rand( 0,1) *(*”“-x™ ln ) (11) 

■'I c=o ] J J 

j =1,2,...,« e , i=\,2,...,N P 

Where rand(0, 1) is a random number between 0 and 1. 

Step 2: Mutation operator: In mutation step, for each 
individual x, (target vector) of the new population, three 
different individuals x r/ , x r2 , and x ri (rl t r21 r3t i ) are 
pseudo-randomly extracted from the population to generate a 
new vector as: 

Zji=x M +F*(x jr 2 -x jr3 ) j = 1,2,..., n c (12) 


ip i +f J P j =P d +PLm (5) 

,=l 7=1 

n b n k 

Y J H j + J j H k =H d ( 6 ) 

7=1 k =l 

/» min < /> < 7] max , i = l,---,n p (7) 

PT" {Hj )<Pj< PT ( H j )• j = !.•••. n b (») 
H™iP ] )<H j <H'r(P J ), j = I,---, n h (9) 

H'r<H k <H'r, k = l,---,n h (10) 

Where P d and H d are system power and heat demands, 
respectively. In general, the transmission losses (PLm) of (5) 
can be calculated through the power generation of all units 
which is known as B-matrix approach (kron’s loss formula 
[1]). P nl,n and P max are unit power capacity limits; h mm and h max 
are unit heat capacity limits. 

The CHPED problem clearly introduces the complication 
of more constraints than in required pure power economic 
dispatch problem. The insufficiencies difficulties with 
conventional methods thus follow from the fact that CHPED 
is a highly constrained optimization problem. 

III. The proposed HDE-MU 

A. Differential evolution (DE) 

The DE algorithm is one of the population-based 
optimization algorithms. The steps for implementing DE are 
as follows [21]: 

Step 1: Initial population: A population of N P initial 
solutions randomly distributed in the n c dimensional search 
space of the optimization problem, are initiated. The DE uses 


Where F e [0,2] is a uniformly distributed random number 
which controls the length of the population exploration vector 

(x r2 - x r3 ). 

Step 3: Crossover operator: After mutation step, the 
crossover operator, according to the following equation, is 
applied on the mutation vector Z /V and the vector x p to generate 
the trial vector U jn for increasing population diversity of the 
mutation vector. 


J Zp, if rand (0,1) < CR 
1 Xj t , otherwise 

j = l,2,...,n c , i = 1,2,. ...Np 


(13) 


Where CR e [0,1] is known as the crossover rate which is a 
constant. 

Step 4: Selection & evaluation operator: The selection & 
evaluation process is repeated for each pair of target/trial 
vectors using the evaluation function F( IJ p ) to compare with 
the evaluation function value F(xji), and the better one will be 
selected to be a member of the DE population generation for 
the next iteration (x^ +1 ). 

B.HDE 

In HDE [22- 26], the one-to-one competition will have a 
faster convergence speed to give a higher probability toward a 
global optimum with much less computation time. It can use a 
small population in the evolutionary process to obtain a global 
solution. Generally, Evolutionary optimization involves two 
critical issues evolutionary direction and population diversity. 
As the evolutionary direction is effective in searching, the 
strong evolutionary direction can reduce the computational 
burden and increase the probability of rapidly finding an 
(possibly local) optimum. As population diversity is increased, 
the genotype of the offspring differs more from the parent. 
Accordingly, a highly diverse population can increase the 
probability of exploring the global optimum and prevent a 
premature convergence to a local optimum. These two 
important factors are here balanced by employing the 
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migration [27] into HDE that can determine an efficacious 
direction in which to search for a solution and simultaneously 
maintain an appropriate diversity for a small population. The 
migrant individuals are generated based on the best 
individual, x° + ' = ,x^fi' ) , by non-uniformly 

random choice. Genes of the i‘ h individual are regenerated 
according to 



where j =1, ..., n c ; i =1, ..., N p , r k and p are random numbers 
in the range of [0,1]. The migration may be performed if only 
the best fitness has not been improved for over 500 
generations running, and the migrant population will not only 
become a set of newly promising solutions but also easily 
escape the local extreme value trap. More details of the HDE 
have shown in [22— 26]. 

C.MU 

Michalewicz et al. [28] surveyed and compared several 
constraint-handling techniques used in evolutionary 
algorithms. Among these techniques, the penalty function 
method is one of the most popularly used to handle 
constraints. In this method, the objective function includes a 
penalty function that is composed of the squared or absolute 
constraint violation terms. Powell [29] noted that classical 
optimization methods include a penalty function have certain 
weaknesses that become most serious when penalty 
parameters are large. More importantly, large penalty 
parameters ill condition the penalty function so that obtaining 
a good solution is difficult. However, if the penalty 
parameters are too small, the constraint violation does not 
contribute a high cost to the penalty function. Accordingly, 
choosing appropriate penalty parameters is not trivial. Herein, 
the MU [30] is introduced to handle this constrained 
optimization problem. Such a technique can overcome the ill 
conditioned property of the objective function. Considering 
the nonlinear problem with general constraints as follows: 
Where h k (x) and g k (x) stand for equality and inequality 



Fig. l The flow chart of the HDE-MU 

* N 0 : maximum number of iterations 
of outer loop 

**Nj : maximum number of iteartions 
of inner loop 


min /(x) 

X 

subject to h k yx) = 0, k = \,...,m e (15) 

g k (x)<0, k = l,...,m, 

constraints, respectively. 

The augmented Lagrange function (ALF) [29] for 
constrained optimization problems is defined as: 

L a (x, v, o) = f (x) + X a k { [h k (x) + v k f - v* } 

k =i (16) 

+ EA-{ (g k to + <4 )+ - Vk } 

k= 1 

where a k and J3 k are the positive penalty parameters, and the 
corresponding Lagrange multipliers v = (v 1 ,...,v m ) and 

u = (i\,...,u m ) > 0 are associated with equality and 

inequality constraints, respectively. 

The contour of the ALF does not change shape between 
generations while constraints are linear. Therefore, the 
contour of the ALF is simply shifted or biased in relation to 
the original objective function, fix). Consequently, small 
penalty parameters can be used in the MU. However, the 
shape of contour of L a is changed by penalty parameters while 
the constraints are nonlinear, demonstrating that large penalty 
parameters still create computational difficulties. Adaptive 
penalty parameters of the MU are employed to alleviate the 
above difficulties. More details of the MU have shown in 
[30], 

D. The proposed HDE-MU 

Figure 1 displays the flow chart of the proposed algorithm, 
which has two iterative loops. The ALF is used to obtain a 
minimum value in the inner loop with the given penalty 
parameters and multipliers, which are then updated in the 
outer loop toward producing an upper limit of L a . When both 
inner and outer iterations become sufficiently large, the ALF 
converges to a saddle-point of the dual problem [28]. 
Advantages of the proposed HDE-MU are that the HDE 
efficiently searches the optimal solution in the economic 
dispatch process and the MU effectively tackles system 
constraints. 

Rower (MW) 



Rower (MW) 



Fig. 2 Heat-power feasible operation region 
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IV. System Simulations 

This section considers a practical CHPED system which 
contains four conventional power units, two cogeneration 
units, and one heat-only unit. The system power demand (Pf) 
and heat demand ( H d ) of this CHPED problem are 600MW 
and 150MWth, respectively. The valve-point effects and 
transmission loss are both considered, and this example used 
the linear boundary of the Heat-Power feasible region for 
cogeneration units. The associated data of this CHPED 
problem are the same as [19]. 

Figure 2 clear shows the feasible power and heat output 
region of cogeneration unit 1 and 2, and demonstrates 
correlations between the power and heat outputs. 

The MU algorithm was used in HDE to hand the equality 
and inequality constraints. The computation was implemented 
on a personal computer (CPU clock 2.0GHz with 4G Ram) in 
FORTRAN-90 language. Setting factors used in this example 
are as follows; the population size N p is set as 5 for HDE-MU. 
The iteration numbers of outer loop and inner loop are set to 
(outer, inner) as (20, 3000). The implementation of the 
proposed algorithm for this example can be described as 
follows: 

L a {x,v,v) = f(x) + £«*.{ [h k (x) + v k f -vl } 

*=1 

+za {(**(*)+«*)!-°*} < 17) 

k =1 

objective : min f(x) 

A=(Pj ,p 2 ,p 3 ,p 4 ,p 5 ,p 6 ,h 5 ,h 6 7 ) 

fix) = X Cost, (P t ) + £ CostJ (Hj ,Pj) + Cost 7 (tf 7 ) (18) 
t=i i =5 

where 


Cost x (P{) = 0.008/f + 2 P x + 25 + |l00sin(0.042 * (10- P t )| 
Cost 2 (P 2 ) = 0.003 P 2 +1.8P, + 60 + |l40sin(0.04 * (20 - P 2 )\ 
Cost 3 (P 3 ) = 0.0012P, 2 + 2.1 P 3 +100 + |l60sin(0.038*(30 - P 3 )] 
Cost A {P A ) = 0.001P 4 2 +2P 4 + 120 + ]l80sin(0.037 *(40-P 4 ))| 

Cost 5 (H 5 ,P 5 ) = 0.0345Rf + 14.5P 5 + 0.03 H 3 + 

4.2H 5 +0.031P 5 H 5 +2650 
Cost 6 (H 6 ,P 6 ) = 0.0435P 6 2 + 36P 6 + 0.027 h\ + 

0.6H 6 + 0.01 1P 6 H 6 + 1250 
Cost 1 (H 1 ) = 0.038H 2 + 2.0109H 7 +950 

Subject to; 

h \: P 1 + P 2 + P 3 + P 4 + P 5 + P b — P A - PLm = 0 

ft 2 : H 5 +H 6 +H 1 ~H i = 0 

g,: 1.78191489361// 5 - P 5 - 105.744680851 < 0 

g 2 \ 0.177777778// 5 + P 5 -247.0 <0 

g 3 : -0.16984732824// 5 - P 5 +98.80 <0 

g A : 1.15841584158// 6 - P 6 - 46.8811881188 < 0 

- P 6 + 44.0 < 0, if H b < 15.9 
- 0.06768189509// 6 - P 6 +45.0761421319 <0, if H b > 15.9 
P b - 125.8 <0, if H b < 32.4 
0.15116279069// 6 + P b - 130.697674418<0, if // 6 >32.4 

This minimum cost problem consists of one objective 
function with nine variable parameters, ( P 2 , P 2 , P <, P4, P5, Pa, 
H 5 , H f „ and H 7 ), two equality constraints, (h/ and h 2 ), and six 
inequality constraints, (gj~ g 6 ), in the Heat-Power available 
region. 

This CHPED example is used to illustrate the effectiveness 
of the proposed HDE-MU with respect to the quality of the 
solution obtained. This test shows the validly and 
practicability of general use for the CHPED problems. Table 
1 compares nine computational results obtained from the 


Table 1: Comparison of the proposed HDE-MU with previous methods 


Methods 

Items 

HDE-MU 

GSA [2] 

OGSO [9] 

GSO [4] 

OBTLBO 

[18] 

TLBO [18] 

TVAC-PSO 

[19] 

KHA [20] 

CPSO [19] 

P( Gl) 

45.56400 

48.7638 

45.6672 

45.6188 

45.8860 

45.2660 

47.3383 

46.3835 

75.0000 

P( G2) 

98.53982 

98.7469 

98.5401 

98.5401 

98.5398 

98.5479 

98.5398 

104.1223 

112.3800 

P(G3) 

112.67349 

112.0000 

112.6735 

112.6727 

112.6741 

112.6786 

112.6735 

64.3729 

30.0000 

P(G4) 

209.81582 

208.5113 

209.8158 

209.8154 

209.8141 

209.8284 

209.81582 

246.1853 

250.0000 

P( G5) 

94.14597 

92.6909 

94.0532 

94.1027 

93.8249 

94.4121 

92.3718 

98.9736 

93.2701 

P( G6) 

40.00000 

40.0000 

40.0000 

40.0001 

40.0002 

40.0062 

40.0000 

40.7401 

40.1585 

H( G5) 

27.40126 

35.9704 

27.9472 

27.6600 

29.2914 

25.8365 

37.8467 

0.0000 

32.5655 

H( G6) 

75.00000 

75.0000 

74.9996 

74.9987 

75.0002 

74.9970 

74.9999 

66.7100 

72.6738 

H( G7) 

47.59874 

39.0000 

47.0532 

47.3413 

45.7084 

49.1666 

37.1532 

83.2900 

44.7606 

SP(G) (MW) 

600.73910 

600.7128 

600.7498 

600.7498 

600.7391 

600.7392 

600.73922 

600.7777 

600.8086 

PLm (MW) 

0.73910 

0.7128 

0.7498 

0.7497 

0.7394 

0.7392 

0.73922 

0.7777 

0.8086 

G) 

(MWth) 

150.00000 

149.9704 

150.0000 

150.0000 

150.0000 

150.0001 

149.9998 

150.0000 

149.9999 

TC ($/h) 

10094.21766 

9912.6928 

10094.2407 

10094.2670 

10094.3529 

10094.8384 

10100.3164 

10111.1501 

10325.3339 

CPU time 
(Sec) 

2.507 

2.578 

2.5681 

2.4203 

3.06 

2.86 

3.25 

N/A 

3.29 
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Fig. 3 Convergence characteristics of the proposed HDE-MU 
in comparison with DE-MU. 


proposed HDE-MU, gravitational search algorithm (GSA) [2], 
opposition-based group search optimization (OGSO) [9], 
group search optimization (GSO) [4], oppositional teaching 
learning based optimization (OTLBO) [18], TLBO [18], 
particle swarm optimization with time varying acceleration 
coefficients (TVAC-PSO) [19], krill herd algorithm (KHA) 
[20] and classic PSO (CPSO) [19]. Although results of table 1 
all located in the Heat-Power feasible operation region, the 
results of GSA [2], TLBO [18], TVAC-PSO [19] and CPSO 
[19] are infeasible solutions, because their heat outputs didn’t 
match the heat demand (H d = 150MWth) that will cause an 
unbalance of system’s heat demand. By investigating the 
results presented in Table 1, it is observed that the best total 



No. of trials 

Fig. 4 Distribution of the minimum costs in 50 trials. 

$/h, which is much less than the best results previously 
reported in TVAC-PSO [19], KHA [20] and CPSO [19], and 
also slightly better than the best results of OGSO [9], GSO [4], 
OBTLBO [18] and TLBO [18]. As seen from the best 
solution of HDE-MU listed in Table 1, the power output is 
600.73910 MW, the heat output is 150 MWth and the 
transmission loss is 0.73910MW. Apparently, the equality 
constraints of heat and power demands are fully satisfied. 
Therefore, the result obtained from the proposed HDE-MU is 
an optimal and feasible solution. According to the absolute 
CPU time of different methods listed in the Table 1, the 
proposed HDE-MU is faster than most of the compared 
algorithms, but slightly slower than GSO [4]. Consequently, 
the HDE-MU is not only supplying the exact total of heat and 
power, but also has the less CPU time than most approaches 


for the cogeneration system. Table 1 reveals that the 
HDE-MU is a satisfactory algorithm for solving the CHPED 
problem. The convergence characteristics of the proposed 
method in comparison with DE with MU (DE-MU) for this 
real example are depicted in Fig. 3, and the distribution of the 
minimum costs after 50 trials by the HDE-MU is plotted in 
Fig. 4. It can be observed that the worst, average and best 
costs are very close. The small standard deviation again 
confirms the stability of the proposed HDE-MU. 

V. Conclusions 

The HDE-MU for solving the CHPED problem has been 
proposed herein. Complication of the CHPED problem lies in 
the constraints imposed by the mutual dependencies of 
multi-demand and Heat-Power capacity. The HDE helps the 
proposed method efficiently search and refined exploit. The 
MU helps the proposed method avoid deforming the ALF and 
resulting in difficulty of solution searching. The proposed 
algorithm integrates the HDE and the MU that has the merits 
of taking a wide range of penalty parameters and a small 
population. A practical CHPED system with transmission 
losses and valve-point effects is used to compare the proposed 
HDE-MU with eight reported methods. Simulation results 
demonstrate that the proposed algorithm has more advantages 
for solving the CHPED problem than the previous methods. 
The contributions of this paper are the MU effectively handles 
the Heat-Power feasible region constraints, the HDE 
efficiently searches the optimal solutions in the economic 
dispatch process, and the author provides a valid and 
efficacious algorithm for the CHPED problem. 
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